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AbptMet: Two elliptic grid generation schemes thnt yield orthogonnl 
grids nre presented (FGBR and UBGR). Each scheme is different from the 
original elliptic grid generation scheme developed by Thompson, et al., in 
the treatment of boundary points^ With the UBGR scheme, the grid points 
on the flow boundaries are automatically determined by the algorithm, 
while with the FBGR scheme at most one half of the boundary grid points 
may be prespecilled and the remainder of boundary grid points are deter- 
mined automatically. Numerical examples show their capability of easy 
stretching, clustering and shock fitting while maintaining orthogonality of 
grid. The present method can be implemented into existing elliptic grid 
generation programs with relatively minor modifications. 


§1. Introduction 

The elliptic grid generation scheme developed by Thompson, et al.[l], 
has been regarded as a method to generate non-orthogonal grids. However, 
approximate orthogonality of grids in the \lclnity of flow boundary may be 
achieved by the algorithm devised by Steger and Sorenson [2]. 

Analytical conformal mapping [3], or its numerical version, generates or- 
thogonal grids. The advantages of completely orthogonal grids are (a) the 
difference equations on the transformed coordinates are simpler than on non- 
orthogonal grids, and thus both computing time and memory space may be 
saved, (b) accuracy of the difference equations on orthogonal grids is better 
than on non-orthogonal girds, (c) if the computational domain is rectangular, 
a fast direct solver may be applicable in place of iterative schemes. 

The most fundamental difference between the conformal mapping and the 
elliptic grid generation method is that the grids on the flow boundary are 
prespecified in the elliptic grid generation, while in the former the grids on 
the boundary are determined py the confomal mapping. 



The objective of the present paper is to show that orthogonal grids may 
be generated by using the elliptic grid generation equations by automatically 
determining at least one half of the boundary grid points as part of the grid 
generation. The present method is more versatile than conformal mapping 
because (a) clustering and stretching of grids are easy, (b) adapting to the 
flow solution is possible, (c) orthogonal grids may be generated even when 
one half of the bo«iudary grid points are flxed, and (d) the algorithm may be 
extended to control orthogonality of grids in three-dimensional elliptic grid 
generation. 

The coefficient of the cross derivative term in the elliptic partial diflTerential 
equations on the computational dom.^ m, which will be called the B-coefficient, 
becomes zero if the transformation is orthogonal. The simpliflcation of the 
flow equations on orthogonal grids is due to this fact. However, an orthogonal 
coordinate transformation does not necessarily mean that the B-co^'tfficient in 
numerical computations is exactly zero [4). This is because even if the coor- 
dinate transformation is analytically orthogonal, the B-coefficient in numeri- 
cal computations is calculated by finite diflTerence approximation. Therefore, 
from the view point of numerical computations^ it is desirable to have the 
numerically calculated B-coefficient exactly zero, than to have analytically 
orthogonal grids as the case of conformal mapping. In this paper we use 
the term numerical orthogonality” meaning that the B-coefficient calculated 
numerically is zero. Although the grids generated by the proposed shemes 
are not numerically orthogonal in this sense, but are rather approximation to 
analytically orthogonal grids. However, feasibility of generating numerically 
orthogonal grids by using elliptic systems is discussed later. 

In the remainder of this paper, (1) the elliptic PDEs and the boundary 
conditions to generate analytically orthogonal grids are derived from or- 
thogonality conditions, (2) the elliptic equations and boundary conditions 
are transformed onto the computational domain, (3) difference equations for 
the elliptic grid generation equations and boundary conditions are derived, 
(4) grids generi^ted by the present method for a number of geometries are il- 
lustrated, and finally (5) feasibility of generating numerically orthogonal grid 
is discussedf 

§2. Elliptic Partial Differential Equation! and Boundary Conditions for 
Orthogonid Coordinate Transformations 

Denoting s-y and as physical and cwputational coordinates, respec- 
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tlvely^ let us consider two families of curves^ Hx,y) = w and ri(x,y) = d, 
on the physical domain as shown in Fig. 1, where c and d are parameters. 
We now determine the conditions for the two families of curves to become 
orthogonal. Along a curve, ( c, we have 

d( = (gdx + iydy = 0 (1) 

or equivalently, 

(dy/dx)i = ( 2 ) 

where subscript 1 indicates that the derivative is the slope of ( = c on the 
x~y plane. Similarly, the slope of the curve tj = d is written as 


(dy/dx)2 = (3) 

The orthogonality between f = c and r) ^ d requires that Eq.(2) be the 
negative reciprocal of Eq.(3), that is, 

^ = -^ = ak(x,y) (4) 

kti f* 

where a is a constant, and k{x, y) is a continuous function which may be 
arbitrarily specified. Equation (4) may be equivaiently written as 

tfx — O'Hy vgs 

The transformation between x-y and is conformal if ak = 1. Otherwise 
the transformation is orthogonal but not conformal. 

Equation (5) is a sufficient condition for orthogonaiity, but it is not suitable 
in determining ^( 3 ;, y) and rj{x, y). To find more appopriate equations, we 
eliminate ( or rj from Eq. (5) (one at a time) and obtain elliptic PDEs as 
follows 


1 1 
i'j^Vx)x “f* ~ ® 
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Notice that the above equations do not involve the constant a. Equation 
(6) is a convenient form b^ause the equations are elliptic PDEr^ for which 
there are many solution methods available. Unfortunately, Eq. (6) is not 
a sufficient condition for orthogonality any more, but is only a necessary 
condition. This is because, while Eq. (6) is derived from Eq. (5), the reverse 
is not possible without additional conditions. However, the solution of Eq. 
(6) are uniquely determined when boundary conditions are specified. In fact, 
the orthogonal set of ( and tj are determined with appropriate boundary 
conditions. 

In order to derive boundary conditions rigorously, we write the following 
theorem: 

[Theorem] Suppose ^{x, y) and ri{x, y) are smooth functions on the. x-y plane. 
In order that the family of curves r)(x,y) = d are orthogonal to a curve 
((x, y) = c, the partial derivative of ff normal to the curve f = c must be 
zero. 

[Proof] The gradient of the curve f = c is given by Eq.(l). The derivative of 
t/(x, y) normal to f = c may be written as 


y)!Sn = + nyh (7) 

where lx and ly are directional cosines of the outward normal to f = c. 
Referring to Fig. 2, the gradient of f = c is related to the directional cosines 
as 


(dyldx)i=^-lxlly (3) 

Eliminating (dy!dx)\ with Eq.(2) and a slight rewriting yields 

^y/ ^y “ (fi) 

Notice that, if ly approaches zero, fy also approaches zero = c becomes 
vertical) while and are both nonzero. Therefore, even when ly appraches 
zero, the left side of Eq.(9) remains nonzero and nonsingular. 

Introducing Eq.(9) into Eq.(7) yields 

dri(x, y)ldn = {lylQinxix + rjyiy) ( 10 ) 
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tf the left side of Eq.(10) is set to zero, then, 


“f* = 0 (11) 

which is equivalent to Eq.(4). Thus ^ c and rf = dm orthogonal. 

[CoroUaryj If a boundary of the ar-y domain is represented by a curve f = r. 
the derivative of ri{x,y) normal to ( === 0 must be zero along that part of 
the boundary. Similarly, if a boundary of the x-y domain is repreiiented by 
a curve rj = d, the derivative of f(x, y) normal Xotfs^d must be zero. 

The above corollary may be intepreted as that, when the value of f (or rj) is 
set to a constant along a portion of the flow boundary, the normal derivative 
of 1) (or ^ respectively) must be zero along that portion. 

In applying the above rule to grid generation, suppose the whole boundary 
of the x-y domain is divided into four parts, Gg, y = 1 to 4, each of which 
is transformed respectively to the left, top, right and bottom sides of the 
rectangular computational domain, 0 < ( < (max ai^d 0 < ri < rtmax (see 
Fig.3a). Thus, the boundary values of f and rj are specjfled as 


( = 0 along Gi 
V “ Vmax uioilg G2 
< “ (max along G3 
77 = 0 along G4 


( 12 ) 


The additional boundary conditions are set in accordance with the corollary 
as follows: 


drf(x,y)/dn = G along Gi and G3 
d((x, y)ldn — 0 along G2 smd G4 

In conclusion, an orthogonal tranformation of the x-y domain to the (-}] 
domain is completed if Eq. (6) with bounday conditions Eqs.(12) and (13) are 
solved. 


§3. Derivation of Grid Generation Equations 

The grid generation equations that map the uniformaly spaced Cartesian 
grids on the computational domain onto the physical domain are obtained 
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by inveriinK Bq.( 6 ) and the boundary condition*, Ev.(12Hl3), in accordance 
rritb 


C» = VnN 
e» = -tnN 

n* = -vt/J 

*= 

where J is the Jacobian given by 

J = *<y, - z„y( 

Equation ( 6 ) is transformed to 

- Vifii-* = ~(J /k)]v^k( - y<y,*, +x^kf- *«*,*,] 
X(Ly — yi(Lx == -(//fc)IyjA;, — y,y«*< + if*,, - x^X(k(] 


(14) 


(IS) 


In the above equations^ L is a quasi-linear elliptic partial differential operator 
given by 


where 


L^A^-2B 


d* a® 

d(dti ^ dt)X 


B= XfX,, + y(y„ 

C = 


Separating Ly and Lx in Eq.(15) yieids 

Lx = ( 07 + 2 ®^®,/?) 


where 


« = (Vfikf - y(k„) 

P = (®e*< - *<**) 

7 = (®ey* + **%) 


(16) 


(17) 
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The boundary conditions^ Eq.(12), are implemented on the computational 
domain as follows: First, define 


/i(a:,p) = 0 on ^ — 1,2,3 and 4 (18) 

where , ^ = 1, 2, 3 and 4, are the left, top, right and bottom boundaries 
of the computational domain respectively, and /g = 0 is the shape of the 
corresponding boundary on the x-y domain. By using Eqs.(7) and and (14), 
Eq.(13)becomes 


= 0 along G,(f = 0) and Gs(f — (ma() 

= 0 along Gi{fj = and G 4 (r) = 0) ^ ' 

where Ig and ly are directional cosines of the outward normal on the boundary. 
A central difference scheme is used to discretize Eq.(17) to yield 


A(xi^ij - 2xi,j + H- - 2ar,j + a?i,,.|.i) 

= B{xi^ij^i — Xi—ij^i — Xi^ij^i + ^»— + 

( 20 ) 

Mvi—lj + j/i-f l,i) + + 2/i,j+l) 

“ F(yt‘-|-ijj-|-i ““ y % — — 1 “4^ yi — i,j — 1)/2 ”1“ 

= ft) 

where R[^j and represent the difference approximations for the right 
sides of Eq.(17). 

Equation (20) is solved iteratively. Any of the iterative methods used for 
existing elliptic grid generation scheme should work for the present method. 


§4. Undetermined Boundary Grid Relaxation (UBGR) icheme 

The numerical treatment of boundary conditions is different in the two 
schemes proposed: UBGR and FBGR. In the UBGR scheme, all the grid 
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points on the boundary are undetermined, while up to one half of the boun- 
dary grid points may 1^ prespecifled in the FBGR scheme and the remaining 
boundary grid points are undetermined. 

Detail of the treatment of the boundar| 1;onditions for the IJBGE scheme 
is described in the remainder of this section. The boundary conditions for 
the FBGR scheme are explained in a later section since it is an extension of 
the UBGR scheme. 

For internal grid points, the two equations in Eq.(20) detemine Xij and 
yi, 3 'f provided that the coordinates of the surrounding eight grid points are 
given. Of course the coordinates of the surrounding grid points are iteratively 
revised in every iteration step. However, assuming temporarily that the 
coordinates of the surrounding grid points are known simplifies explanation 
of the derivation of the diflSerence equations for the boundary grid points. 

For a grid point on the boundary, the two equations in Eq.(20) can not 
be satisfied simultanously because of the following additional constraints: (1) 
boundary grid must be on the specified boundary line, and (2) orthogonality of 
grid lines at the bounday is required. Althogb the othogonality requirement 
can be Incorporated into the elliptic equation as a bounday condition, the 
first constraint does not allow the two elliptic equations to be independent. 
This means that only one of the two equations can be used unless the two 
equations are linearly dependent. 

In order to find an appropriate equation, consider a boundary grid point 
{i,j] to be determined on a boundary as depicted in Fig, 4a. We assume that 
the three surrounding grid points denoted by L, R, and B are temporarily 
fixed. The variables s and t are local Cartesian coordinates as shown in Fig. 
4b. The elliptic equations for x and y may be written in terms of s and t as 


Asii - 2BS(„ + Cs„ = 

(21a) 

+ Ct,, = hC) 

(2Ib) 


Since Eqs.(21a) and (21b) are independent, one has to choose only one of 
the two equations. We choose Eq.(21a) and abandon (21b) for the following 
reasons: The coordinate t is in the normal direction froin the boundary so ti^ j 
on the boundary is not allowed to change. On the other hand, s is tangent 
to the boundary so has a freedom to change. 


8 


The difference form of the flrtt and third terms of Eq.(21a) may be written 
as 


•u = - 2«<,> + 

(22a) 

“ (*i>j ■" i) 

■" 1/2 

(22b) 


The orthogonality boundary condition requires that the first term on the right 
side of Eq.(22b) be zero, 


® ( 23 ) 

By introducing Eq.(22) into Eq,(21a) and using Eq.(23), we obtain 

1 j) + (24) 

where B in Eq.(21a) is set to zero because of orthogonality at the boundary. 
If the coordinates of the surrounding grid points are known, Sij is found by 
solving the above equation. Then U^j is calculated by introducing the value 
of Sij into the equation that defines the boundary. 


Mf) = o 


(25) 


§5r Control of Grid Spacing 

Control of grid spacing for clustering and stretching is achieved by the 
following two algorithms: (1) variable grid spacing on the computational grids, 
and (2) space dependent coefficient^ k{x, y). 

In the first approach, the Cartesian grid points with variable spacing on the 
computational domain are mapped onto the physical domain as illustrated in 
Fig. 5. Small grid spacing on a part of the physical domain is obtained by 
imposing a small grid spacing on the computational grid in the corresponding 
area. This technique may be used on the ( and rf directions independently. 
Once the coordinates of the grid points on the physical domain are deter- 
mined, the grid may used as if it were generated for an uniformely spaced 
computational grid. 
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The second approach is to change k In space, Since k in Eq.(4) is propo- 
tional to the aspect ratio of a grid cell, the clustering effect in the i^-direction, 
for example, Is obtained by reducing the value of k In the interested area 
on the physical domain. While the first approach can only change the grid 
spacing uniformly In the ( or ^ direction, the second approach can cluster 
or stretch grid more fiexibly on the physical domain. Uses of k{Xj y) are 
illustrated later with numerical examples. 

In addition, these two techniques may be applied simultaneously in one 
grid generation as needed. 

§6. Fixed Boundary Grid Relaxation (FBGR) Scheme 

In the TJBGR scheme two end points of a grid line are undetermined so both 
move along the physical boundary during iterations, while the grid line on 
the computational domain is fixed. Amother way of satisfying the orthogonal 
grid generation equations is that: (a) one end point of a grid line is fixed, (b) 
the other end point of the grid line is left to be variable, and (c) grid spacings 
adjacent to the corresponding grid line on the computational domain are 
variable. This latter approach is denoted as the FBGR scheme. 

This scheme can be implemented by modifying the UBGR scheme only in 
the treatment of the fixed boundary grid points as explained next. 

Suppose that the grid points with black circles in Fig. 6 are obtained after 
an iteration step with a given set of grid spacings 6(i s (i* — and 
— (i. First, the hypothetical coodinates to satisfy Eq. (24) 
are calculated (open circle marked as &) as if the grid point ( 1 , j) on the 
boundary is an undetermined point. Second, the following two chord length 
ratios are calculated: 


r AB/(AB + BC) ^ ^ 

where ABf BC^ ABf and &C are chord lengths. Third, the computational 
grid spacings and are replaced by 



(2T) 


(28) 
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where 


a = tf e./r' - ife+iAr* - 1) 
and 


6 = { 6 (i + 6 (i^i) — a 

The geometrical interpretations of the above equations are : 

(i) If and 6 (i^i are the computational grid spacings, the orthogonal 
boundary condition requires that the boundary grid must be at 

(ii) If the boundary grid is to be located at B while satisfying the or- 
thogonality boundary conditions, the computational grid line i must move 
toward liiio i — 1. 

(iii) The functional relation among 6 (i^i and r' are fitted by a quadratic 
polynomial, and then the revised value of the grid spacing 6(4 is calculated 
by the quadratic function thus obtained. 

(iv) The total computational grid spacing between (1 — 1) and (1 + 1) is 
unchanged through this modification, namely 

+ + ( 29 ) 

grid spadngs next to all the fixed boundary grid points are modified 
after each iteration step. 


§7. Numerical niuftratioDi 
Example 1: 

Figure 7 shows the grid generated by the preisent method for a quadrilateral 
domain. An uniformly spaced Cartesian grid on the computational domain 
is mapped to the the physical domain by the UBGR method with k ^ I, 

Example 2 : 

Figure 8 shows the grid generated for a NACA-0012 airfoil by the UBGR 
scheme with k and no clustering. The grid is not, however, orthogonal 
in the vicinity of the trailing edge since the two grid points at the trailing 
edge are fixed without any FBGR treatment. 

Example 3: 
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Figure 9 shows the grid generated for & NACA-0012 airfoil by the UBGR 
schem with a clustering effect created by a variable grid spacing on the 
intermediate computational domain. 

Eiample 4: 

Figure 10 shows the grid points clustered around a prescribed curve to 
illustrate the capability of solution adapting. This was done by setting k{x, y) 
in Eq.(4), or equivalently £q.(18), as 


where ci and C 2 are positive constants and h{x, y) = 0 is the curve (for 
example, a shocjk profile) to which grid points are to be clustered. The 
function h{x, y) used in the present example is 


h(Xf y) sinix) — y 


Example 5: 

Figure 11 s’aows the C-type grid generated by the UBGR scheme for a 
cambered air/oil. The grid points on the Kutta cut do not match because the 
points above and below were determined independently. 

Example 6; 

Figure 12 shows the C-type grid generated by the FBGR for NACA-0012 
airfoil. The grid points on the airfoil surface and the Kutta cut are fixed, 
while the grid points on the outer boundaries are all determined by the UBGR 
scheme. The grid points near the trailing edge are clustered for this particular 
grid but this clustering is not essential to the present grid generation scheme. 

Example 7: 

Figure 13 is for the same configuration as Example 6 except the airfoil is 
cambered. 

All the grids shown above were calculated by the point successive over 
relaxation iterative scheme, which was used for simplicity of programming. 
Although no effort has been paid in the present study to accelerating the con- 
vergence rate, the convergence of the iterative scheme is felt to be slow par'' 
ticularly in comparison with the conventional elliptic grid generation scheme 
with fixed boundary grid points. It seems possible to increase the computar 
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tiona) efficieBcy for the proposed scheme by the ADI or multigrid method or 
both. 


§8. NumericaUy Orthogonal Grids 

The orthogonal grids generated by the UBGR and VBGR schemes are 
approximation to analytically orthogonal grids. This is because the basic 
equations are based on analytically orthogonal transfomations but the finite 
difTerence approximation was used to generate grids, In this section the 
possibility of generating numerically orthogonal grids is analyzed. 

The condition of numerical orthogonality may be defined as 


(«jaK*ja:) + (^<yK%) = 0 (30) 

where 64 and 6 j are central difference operators in the i and j directions 
respectively. Rewriting Eq.(30) yields 


diuldjx = -‘dix/bjy = ak(xt y) 
which may be further rewritten as 


6iy = akSjX 
6iX = —akdjp 


By rewriting, we obtain 


(31) 


(32) 


«i(r«<®) + = 0 

K 

H\Siy) + a%ikS,y) = 0 


(33) 


The coefficient a in Eq.(33) is eliminated as follows. The two equations In 
Eq.(32) are added after taking the square of each term, and the resulting 
equation is solved for as 

fl2«C/(A:U) (34) 
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where 


C = (Sixf + (Sivf 


By introducing B<i.(34), Eq,(33) becomes 


Ah,%ilSix) + cls0Six) = O 


(36) 


If the difference operator is replaced by the differential operator^ Eq.(36) 
becomes equivalent to Eq.(17) except Eq.(36) is missing the B-coefflcient* 
This last point is understandable because Eq.(36) has been derived from the 
orthogonality condition. 

Difference equations are meaningless unless the meaning of the difference 
operator 6 is clearly defined. So we define ^ here as the central difference 
operator associated with one grid interval. Then Eq.(36) is in the well known 
five point difference form of elliptic partial differential equations. The five 
points involved are shown with black dots in Fig. 14, which are called here 
the primary grid points. 

With the present definition of 6. however, Eq.(30) may be written more 
explicitly as 


-*(3/i-f-l/2,y+l/2 - y»-.l/2,j+l/2X2/»,i+X Vij) = 0 


(3Ta) 


and 

(«I-H J — Xi,j){Xi^i/2,j^i/2 — Xi^i/2,j^i/2) 

~ 3/i,i)(y»+i/2,j+i/2 — y*+i/2,j-i/2) = o 

where (i -j- 1/2, j + 1/2) is a grid point in the staggered position as shown 
with an open circle in Fig. 15. The grid points at the staggered positions 
are called the auxiliary grid. A system of elliptic difference equations in the 
same form as Eq.(36) can be derived for the auxiliary grid from Eq.(37). The 
important thing to be observed here is that Eq.(37) is interpreted as the con- 
dition for numerical orthogonality between the primary and auxiliary grids, 
but not for the numeical orthogonality among the primary grid points or 
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among the auxiliary gid points.. 5pt Another way of defining the meaning 

of 6 is possible. For example, suppose 6 is the central difference operator 
associated with two grid intervals. Then, the left side of Eq.(30) l^omes ex- 
actly the central difference approximation of the B-coefflcint explained earlier. 
Therefore, based on this definition the grid generated whould be numerically 
orthogonal. The problem, however, is that the same five difference equations 
that generates the primary grid can not be derived any more. 

The present observation leads us to the conclusion that, if the grid is 
generated by the conventional five point difference equation, the B^coefficient 
calculated numerically by the central difference approximation with two grid 
intervals can not be made exactly zero. 

On the orther hand, there Is a possibility that numerically orthogonal 
staggered grid may be obtained. This will require generating both the primary 
and auxiliary grids by two sets of elliptic grid generation equations coupled 
with appropriate boundary conditions. The detail of such a scheme needs 
further research, however. 

§9. Ooneiuiioni 

The present paper showed that orthogonal grids may be generated by 
solving a finite difference approximation for the elliptic partial differential 
equations with either of the UBGR or FBGR schemes. Ihese schemes are 
more versatile than conformal mapping because (1) the grid cell aspect ratio 
may be controlled through the function k(x, y), and (2) they can hi used for 
a local control of orthogonality of grid, and (3) the schemes can be extended 
to three-dimensional grid generation to control orthogonality of grids in one 
or two directions. 

The proposed methods can be implemented in existing grid generation 
programs which are based on the elliptic grid generation equations with 
relatively minor effort. 

One drawback of the proposed shemes is that the iterative convergence 
rate is slower than the conventional elliptic grid generation scheme with fixed 
boundary grid points. The first reason for this is that the boundary conditions 
are Neuman lype in the proposed schemes while the boundary conditions of 
the conventional elliptic grid generation schemes are all Dirichlet type. The 
second reason is that the boundary relaxation schemes are nolinear. Since 
the convergence rate of the proposed schemes ! . strongly dependent of the 
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number of grid points, the multigrid method would be useful to increase the 
computing efficiency. 

As an add!.^ional product of the present worlc, it was shown that the grids 
generated by the flve*point central difference approximation to the elliptic 
PDEs can not be made numerically orthogonal if the B-cefflcient is calculated 
by the central difference apprximation associated with two grid intervals, but 
grids may be made numerically orthogonal between two grid systems that are 
in the staggered position to each other. 
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Figure 7 Orthogonal grid generated for a quadrilateral 
domidn with no grid spacing control 



Figure 8 0-grid for NACA-0012 airfoil 
generated bj the UBGR scheme 
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Figure 9 Dlustration of clustering by variable grid spacing 
on the intermediate compuational domain 
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Figure 10 


Grid on a rectangular domain with adapting 
to a given curve 
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(B) 

Figure 11 (continued) Close view of the leading edge 
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Figure 12 C-grid generated for NACA-0012 airfoil by the PBGR scheme 
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Figure 13 C-grid generated for an arbitrarily cambered airfoil by 
the FBGR scheme 
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Figure 14 Primary grid 
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Figure 15 Auxiliary grid 


30 



